Chapter 3 Data Processing
The goal is to compare rainfall and discharge data in order to compute the following metrics: 1) total quantity of response flow, 2) magnitude of peak response discharge, 3) duration of response, 4) lag to peak time (start time of storm to time of peak discharge), and 5) runoff ratio (depth of Q (mm) divided by the depth of rainfall (mm) for the rain storm associated with the runoff response).
3.1 Precipitation
Before starting this analysis, I identified discrete rainfall events from the tipping bucket rain gauge output using the USDA Rainfall Intensity Summarization Tool (RIST; (ARS 2013)). RIST calculated the following metrics for each event: storm duration (hrs), rain depth (mm); maximum intensity (mm h−1) over 5-, 15-, 30-, and 60-min intervals; and erosivity. Because the stage data is recorded every 15 minutes, I’m using the start time of the maximum intensity 15 minute interval as the storm start time. The function below formats the RIST output into the form I want.
#function to read in rain files
rain_reader <- function(rain_file){
df <-read.delim(rain_file, sep = "", header = TRUE) %>%
select(-c(1,3:7))
}
#function to process the storm output
rain_format <- function(rain_file){
df<- rain_file %>%
select(Max_15_start_date, Max_15_start_time, Precip.in.) %>%
rename(date =1, time=2, p_in =3) %>%
mutate(datetime = paste(date, time, sep = " ")) %>%
filter(!(p_in ==0)) %>%
mutate(p_mm = p_in*25.4) %>%
mutate(DateTime = mdy_hms(datetime, tz = "GMT"))
}
#read in and format michigan ditch and bl4 rain data
mich_rain <- rain_reader("raw_sitedata/michiganditch/michigan_storm.txt")
bl4_rain <- rain_reader("raw_sitedata/bl4/bl4_storm.csv")
mich_storm <- rain_format(mich_rain)
bl4_storm <- rain_format(bl4_rain)3.2 Stream Stage
Michigan Ditch
Michigan Ditch stream stage was recorded by a pressure transducer that doesn’t take into account atmospheric pressure, so the data needs to be corrected using barometic data. The continuous stage data is then calibrated according to the manual stage measurements.
#read in pt data
michiganditch_stage_pt=read_csv('raw_sitedata/michiganditch/michiganditch_stage_composite_pt.csv', show_col_types = FALSE) %>%
rename(Pw_kPa = 3)
michiganditch_stage_pt$DateTime = as_datetime(michiganditch_stage_pt$datetime,tz="GMT", format="%m/%d/%Y %H:%M")
#read in baro data
michiganditch_baro=read_csv('raw_sitedata/michiganditch/michiganditch_baro_composite.csv', show_col_types = FALSE) %>%
filter(!is.na(Pressure_kPa)) %>%
rename(Pa_kPa = 2)
michiganditch_baro$DateTime = as_datetime(michiganditch_baro$datetime,tz="GMT", format="%m/%d/%Y %H:%M")
#join pt and baro and convert to stage in cm
mich_stage = full_join(michiganditch_baro, michiganditch_stage_pt, by="DateTime",all=TRUE) %>%
mutate(stage_kPa = Pw_kPa - Pa_kPa) %>%
mutate(stage_cm = (stage_kPa*101.97162129779)/10) %>%
select(DateTime, stage_cm)
#read in file with manual stage measurements
discharge_michiganditch=read_csv('raw_sitedata/michiganditch/discharge_michiganditch.csv', show_col_types = FALSE)
discharge_michiganditch$DateTime = as.POSIXct(discharge_michiganditch$datetime,tz="GMT",format="%m/%d/%Y %H")
#merge the manual discharge measurements with the sensor time series
michditch_stage=full_join(mich_stage,discharge_michiganditch,by='DateTime') %>%
select(1,2,7,9)
#correct stage based on manual stage measurements
michditch_stage<-mutate(michditch_stage, stage_corr_cm = stage_cm)
#offset adjust
michditch_stage$stage_corr_cm=ifelse(michditch_stage$DateTime >ymd_hms('2021-06-21 10:30:00'),
michditch_stage$stage_corr_cm-1,michditch_stage$stage_corr_cm)
michditch_stage$stage_corr_cm=ifelse(michditch_stage$DateTime>ymd_hms('2021-08-05 11:00:00'),
michditch_stage$stage_corr_cm-2.3,michditch_stage$stage_corr_cm)
#delete sensor spikes
michditch_stage <- michditch_stage %>%
filter(!(DateTime < ymd_hms("2021-01-10 09:30:00") |
DateTime == ymd_hms("2021-06-21 09:45:00") |
stage_corr_cm <= 0))
#plot the sensor and manual stage time series
michiganditch=ggplot()+
geom_line(data=michditch_stage,aes(x=DateTime,y=stage_corr_cm))+
geom_point(data=michditch_stage, aes(x=DateTime,y=manual_stage_cm),color='red')+
labs(x= "Time", y= "Stage (cm)") +
theme_bw()
ggplotly(michiganditch)Figure 3.1: Michigan Ditch stream stage. The red dots indicate manual stage measurements.
Blue Lake 4
The Blue Lake data was collected with a capacitance rod, so doesn’t need to be corrected for atmospheric pressure.
#read in file with manual stage measurements
discharge_bl4=read_csv('raw_sitedata/bl4/discharge_bl4.csv', show_col_types = FALSE)
discharge_bl4$DateTime = as.POSIXct(discharge_bl4$datetime,tzone="US/Mountain",format="%m/%d/%Y %H")
discharge_bl4[1:2,10] = as.POSIXct('2021-06-01 15:30:00')
#read in sensor stage data and format
bl4_stage=read_csv('raw_sitedata/bl4/bl4_stage_composite.csv', show_col_types = FALSE)
bl4_stage$Date = as_date(bl4_stage$datetime, format="%m/%d/%Y %H:%M")
bl4_stage$DateTime = round_date(as_datetime(bl4_stage$datetime, format="%m/%d/%Y %H:%M"), unit = "15 mins")
bl4_stage$DateTime = force_tz(bl4_stage$DateTime, tzone = "US/Mountain" )
bl4_stage <-rename(bl4_stage, Stage_mm = wtrhgt__3)
bl4_stage <- bl4_stage[, -c(2:3)]
#convert stage from sensor units to cm
bl4_stage$Stage_cm=bl4_stage$Stage_mm/10
#merge the manual discharge measurements with the sensor time series
bl4_stage=full_join(bl4_stage,discharge_bl4,by='DateTime')
#delete bad measurements-- sensor spikes and other obvious errors
bl4_stage<-filter(bl4_stage, datetime.x != "8/6/2021 8:45")
bl4_stage<-filter(bl4_stage, datetime.x != "6/1/2021 10:33")
bl4_stage<-filter(bl4_stage, datetime.x != "6/1/2021 11:33")
bl4_stage<-filter(bl4_stage, datetime.x != "6/1/2021 12:33")
bl4_stage<-filter(bl4_stage, datetime.x != "6/1/2021 13:33")
#goal is to offset-adjust sensor stage to levels of manual stage
bl4_stage$Stage_corr_cm = bl4_stage$Stage_cm
bl4_stage$Stage_corr_cm = ifelse(bl4_stage$DateTime < '2021-06-01 14:30:00',
bl4_stage$Stage_corr_cm +6, bl4_stage$Stage_corr_cm+2)
bl4_stage$Stage_corr_cm = ifelse(bl4_stage$DateTime > '2021-06-11 12:30:00',
bl4_stage$Stage_corr_cm +4.4, bl4_stage$Stage_corr_cm)
bl4_stage$Stage_corr_cm = ifelse(bl4_stage$DateTime > '2021-07-02 13:45:00',
bl4_stage$Stage_corr_cm +0.2, bl4_stage$Stage_corr_cm)
bl4_stage$Stage_corr_cm = ifelse(bl4_stage$DateTime > '2021-08-06 08:30:00',
bl4_stage$Stage_corr_cm +0.5, bl4_stage$Stage_corr_cm)
#plot the sensor and manual stage time series
bl4=ggplot()+geom_line(data=bl4_stage,aes(x=DateTime,y=Stage_corr_cm))+
geom_point(data=bl4_stage,aes(x=DateTime,y=manual_stage_cm),color='red')+
labs(x= "Time", y= "Stage (cm)") +
theme_bw()
ggplotly(bl4)Figure 3.2: Bl4 stream stage. The red dots indicate manual stage measurements.
3.3 Discharge
Rating curves were developed for both sites using JMP statistical software, and then used to calculate continuous discharge.
#michigan ditch rc is a sqrt transform of the y axis
michditch_stage$Discharge_Ls = (-0.218979 + 0.3533779*(michditch_stage$stage_corr_cm))^2
#plot sensor discharge and manual discharge to make sure rating curve was applied correctly
michiganditch=ggplot(michditch_stage)+
geom_line(aes(x=DateTime,y=Discharge_Ls))+
geom_point(aes(x=DateTime,y=discharge_Ls),color='red')+
xlab("Time")+ ylab("Discharge (L/s)")+
theme_bw()
#bl4 rc is a gompertz 3p sigmoid curve
a=159.07732
b=0.0997502
c=23.905566
bl4_stage$Discharge_Ls = a*exp(-exp(-b*(bl4_stage$Stage_corr_cm - c)))
#plot sensor discharge and manual discharge to make sure rating curve was applied correctly
bl4=ggplot()+geom_line(data=bl4_stage,aes(x=DateTime,y=Discharge_Ls))+
geom_point(data=bl4_stage,aes(x=DateTime,y=discharge_Ls),color='red')+
xlab("Time")+ ylab("Discharge (L/s)")+
theme_bw()
Figure 3.3: Plots of discharge; Michigan Ditch on the left and Blue Lake 4 on the right. The red dots indicate manual discharge measurements.